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What is STECKMAP? 



STECKMAP stands for STEllar Content and Kinematics via Maximum A Pos- 
teriori likelihood. It is a tool for interpreting galaxy spectra in terms of their 
stellar populations, through the derivation of their star formation history, age- 
metallicity relation, kinematics and extinction. From these, a number of in- 
tegrated quantities are also computed, such as luminosity-weighted age and 
metallicity, and mass-weighted age and metallicity (respectively, LWAge, LWZ, 
MWAge, MWZ). 

To do so, the observed spectrum is projected onto a temporal sequence of 
models of single stellar populations, so as to determine a linear combination of 
these models, that fit the observed spectrum best (via a penalized x^). The 
weights of the various components of this linear combination indicate the stellar 
content of the population. This procedure is regularized using various penalizing 
functions. The principles of the method are detailed in [7ll^. 

The STECKMAP software package is public and freely available at 
|http://astro. u-strasbg.fr/~ocvirk/, A number of authors have already 
adopted it and use it in their daily research. Examples of a variety of applica- 
tfons of STECKMAP can be found in H [Til 1 E [10] 

This user guide aims at accompanying the user through the setup and first 
rims of STECKMAP. The last chapter will help the user to understand and 
improve his results and experience with the code. 
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Chapter 1 

STECKMAP setup 



One needs to install Yorick and STECKMAP separately. 

1.1 Installing yorick 

You will thus need to first install Yorick from 

[http://www.maiimae.net/yorick/doc/. This should in principle be painless, 
especially using the pre-compiled binary distributions. Yorick is very light (much 
lighter than the STECKMAP package itself) so you wont have to worry about 
disk space. I recommend installing Yorick in a directory such as $HOME/Yorick/ 
Successful installation should result in the yorick executable to have a full 
path similar to: 

$HOME/Yorick/yorick-2.1/yorick/yorick (for yorick-2.1 that is). This note 
provides few guidelines for installing yorick. Refer to yorick's README for 
more, as it is more complete. On Macs, yorick can be installed through fink. I 
recommend using yorick through emacs. This allows for 2 very cool things to 
happen: 

• enables emacs code editing features: colors, indentation, etc... i.e. a 
helpful coding environment. 

• keyboard shortcuts: the buffer you are currently editing can be executed 
with Ctrl-X Ctrl-S. 

Follow the yorick. el instructions to set this up. You will need to make a 
straightforward modification of your /.emacs. and inside yorick. el dont forget 
to provide the full path to the yorick executable (search for the string "yorick- 
executable-name" and replace the adjacent string "yorick" with the full path to 
the executable). You can check if yorick is correctly installed by launching it 
typing 

bash-3.2$ yorick 

Copyright (c) 2005. The Regents of the University of California. 
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All rights reserved. Yorick 2.1.05 ready. For help type 'help' 

> include, "demo2.i" 

> demo2 

A cool animation should ensue, greeting you to a successful install of Yorick. 



1.2 Installing STECKMAP 

You can find the latest package from here: 

|http : //www . aip . de/Peopre/POcvirk/ Download the latest version. The ver- 
sions are tagged as: 
STECKMAP_laptop_ddmmyy.tar 
or 

STECKMAP_laptop_ddmmyyyy.tar 

In case of problems with the download email me or have a look at the STECKMAP 
webpage: 

|http : //astro . u-strasbg . f r/~ocvirk/ In your Yorick directory (SHOME/Yorick/), 
untar the STECKMAP archive 

tar -xvf STECKMAP. tar 

In order to be able to relocate the code easily, an environment variable 
is used to locate the directory where STECKMAP is installed. Its name is 
STECKMAPROOTDIR. I usually install Yorick and STECKMAP in $HOME/Yorick/, 
so that I can set (for a bash shell): 

export STECKMAPROOTDIR=$HOME/Yorick/ 

In this setup, this rootdir looks like this: 

bash-3.2$ cd $STECKMAPROOTDIR 
bash-3.2$ Is 
STECKMAP 
yorick-2 . 1 



1.3 Checking the successful installation with an 
example 

Launch yorick by typing "yorick" on the shell command line. Once yorick is 
launched, load STECKMAP by typing 

> include, "STECKMAP/Pierre/PDP/sf it . i" 

You will find a couple of example data in Yorick/Pierre/POP/EXAMPLES/. 
They are provided in order to give the user a taste of what can be done and 
how to proceed. Follow this small tutorial. 

STECKMAP is a tool package aiming at constraining the stellar content and 
kinematics from galaxy spectra. To do so, you need two things: 
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• DATA 



• Models against which to compare your data 

The detailed use of the function will be explained in the further chapters. 
For the moment we just want to check STECKMAP is installed properly. 

Data 

See the path and name of an example file I set up by typing 

> fV 

fV is just a variable storing the path and name of the example file. To make 
things easier, it is loaded by default. Convert the example file by typing: 

> a=convert_all(fV) 

this will create the .pdb file, plot the data, and write the redshdift of the 
example galaxy if it is supplied in the fits keywords. Otherwise it will just 
assume 0, as in this case. Note that it will plot the 2d plate and the spectrum 
obtained by summing the whole 2d spectrum in the spatial direction. This is 
the default behaviour of convert_all. 

A model basis is a structure as defined in "Pierre/POP/base^truct.i" . This 
structure contains the sequence of SSPs in time and metallicity. As such, it is 
a data cube. The structure contains this cube and also the wavelengths, the 
resolution, the metallicity scale, and the ages of wach element of the basis. 

SSP basis 

To generate a basis, use the function bRbasisS: 

> b=bRbasis3() 

This will call bRbasisS with all the defaults arguments. To get help and see 
information about the various arguments and their default values you can type 

> help, bRbasisS 
or 

> info, bRbasisS 

Here the default is a PEGASE-HR sequence of SSPs with 10 age bins from 
lOMyr to 20 Gyr with Salpeter IMF and Padova tracks. It is flux-normalizcid by 
default (see the STECMAP paper for details). PEGASE-HR has a resolution 
R=10000 and the wavelength coverage is 4000 — 6900 A. All these values can of 
course be changed when calling bRbasisS. 

To have an idea of the looks of the basis you just generated, type 

> ws 
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this is to clear the display 

> plb, b. f lux (, ,1) ,b. wave 

This will result in a nice plot of the basis for the first metallicity of the basis. 
It can be printed by typing: 

> b.met(l) 

since the mctallicitics arc stored in b.mct. Note; that for computational reasons, 
the metallicities are renormalized. Hence, you can read the original metallicity 
of the first constant-metallicity slice of the SSP cube b.flux(,,l) by typing: 

> Zrescaleml (b .met (1) ) 
it should be 0.05. 

Fitting engine 

The fitting engine is called with the sfit command, which takes as arguments, 
the data, the basis and a bunch of options. Type 

> help, sfit 

to see the current documcmtation for sfit. To fit the data, the structure of which 
is stored in the variable a, using the SSP basis, the structure of which is stored 
in the variable b, type: 

> x=sf it (a,b,kin=l ,epar=3,noskip=l ,sav=l) 

The various options will be explained in the next chapter. 
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Chapter 2 

Fitting your data 



2.1 Preparing your data: convert _all 

Spectra of galaxies come in a variety of forms, ID (such as SDSS), 2D (typical 
long-slit spectroscopy) , 3D (data cubes) , some of them have attached errors and 
masks (e.g. SDSS), others not. STECKMAP on the other hand wiU fit only 
one spectrum at a time. Therefore the first step in the analysis is to convert the 
available data into a ID spectrum in a format that STECKMAP understands. 
This can be done in the majority of cases with the convert_all routine. Since 
even in the case of 3D spectroscopy spectra are collected on a 2D detector and 
even stored in 2D structures (wavelength, fiber number), convert_all can deal 
at the moment only with ID and 2D provided spectra. The inputs and options 
can be obtained as 

> inf o,convert_all 

f una convert _all (f ilelist , cut= , noplot= , log= , zO= , SMRO= , wav= , 
wavaxis= , xs= , xe= , hdu= , f sigm= , errorf ile=) 

A typical call looks as follows: 

> a=convert_all (f V , xs=100 , xe=300 , z0=500 . /300000 . ) ; 
NAXIS2>NAXIS1, assuming AXIS2 wavelength axis 

no redshift in fits header... will assume or take user input 
no S/N in fits header... will assume 100 or take user input 
redshift 0.00166667 

> 

2.1.1 Input 

The only mandatory input is the filename of the data file to be ingested f ilelist. 
This must be a fits file following Note that it can be a list but it is usually more 
practical to convert and analyse spectra one by one. 



8 



2.1.2 Options 



hdu= 

cut= 
log= 

noplot= 



zO= 
SNRO= 
errorf i 



le= 



wavcLxis= 

xs=, xe= 
sigm= 



if fits file contains multiple header data units, specify num- 
ber of hdu to read. 

filelist is cut at cut-th file, default is 10 
log=l enforces log wave sampling in case fits header infor- 
mation is inaccurate. 

disables the plotting of the spectrum read (useful when re- 
motely running a batch of spectra on a machine without 
an active Xll window with nohup for instance, default is 
noplot=0, so plotting happens. 

if rcdshift is not provided in fits header, can be given by 
user as zO= 

same as zO for global signal to noise ratio. Note that ideally 
it is better to supply a noise spectrum via errorf ile=. 
specify a name for an error file in order to fill the sigm 
vector. This will only work for ID spectra right now. 
possible values are 1 or 2. Useful if data is provided as 
2d frame, typical for long-slit spectroscopy. If no value is 
provided (default) we take as the wavelength axis as the 
axis with largest dimension. 

start and end of the stacking in the spatial direction, useful 

for long slit spectroscopy. 

sigm=l forces sigma=l instead of d/(S/N). 



2.1.3 Output 

The direct result of the call is a structure containing various informations about 
the spectrum and the galaxy, and the various filenames of the results files .res, 
the original data file .fits, the intermediate data file .pdb. This latter file is 
written by the routine and will be accessed by the fitting routine sf it. You can 
see the content of the result structure by typing: 

> a 

[galStruct (name="NGC4621" ,f ilename= 

"/Users/Pedro/Yorick/STECKMAP/Pierre/P0P/EXAMPLES/NGC4621-l.pdb",result_dir= 
"Users/Pedro/Yorick/STECKMAP/Pierre/POP/EXAMPLES/",resfile= 

'7Users/Pedro/Yorick/STECKMAP/Pierre/P0P/EXAMPLES/NGC4621-l.res",wavesampling= 
"LIN " , redshif t=0 . 00 166667 , SNR=100) ] 



2.2 Preparing a sequence of SSP models: 
bRbasisS 

There are several versions of the routine for generating the SSP basis. A basis 
is a datacube, the 3 dimensions being (A, age, Z). Note that here in all of 
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Code Reference age range (yr) [Z/H] 

' 'BC03" Bruzual & Chariot 2003 10^ - IJIO^" [0.3,-2] 

"MILES" Vazdekis 2010 2.10^ - l.TlO^o [0.2,-1.3] 

"PHR" Leborgne et al. 2004 2.10^-1.7101'' [0.2,-2] 

"GD05" Gonzalez-Delgado et al. 2005 2.10^^ - 1.710^*' [0.2,-2] 

Table 2.1: References of the SSP models available with STECKMAR 



this user guide the term "basis" is used somehow abusively: because of the 
correlations between the model spectra at different ages and metallicities, such 
a family of vectors is usually not linearly independent and therefore it is not 
a basis in the linear algebra sense. With this warning, we will stick to this 
nomenclature for practical reasons, and because it captures well the idea of 
"projecting" the observed spectrum on a set of model spectra. Note that the 
decomposition obtained in this way is unique only if the set of model spectra 
is linearly independent, i.e. is a basis. The most advanced at the moment is 
bRbasisS. Information can be obtained by typing: 

> info , bRbasisS 

f unc bRbasisS (ages ,&FWHM,&bab,mets=,nbins=,R=,wavel=,N=,basisf ile=,dlEimbda=, 

zr=,br=,inte=,navg=,list=,intages=,dirsfr=,f ilters=,s=) 

> 

A typical call looks as follows: 

> b=bRbasisS([l.e8,2.elO] ,wavel= [S500 . ,7000.] ,basisf ile="MILES" ,nbins=SO) ; 



2.2.1 Inputs 

The most important input here is the age ranges where the models will be 
computed. It is given as a 2-vector [agemin, cigemax]- The unit is years. 



2.2.2 Options 

nbins= Number of ages bins of the basis. 

wavel= Wavelength range of the models (by default the broadest 

available range is taken. 

R= The models can be broadened to an arbitrary spectral res- 

olution to account for instrumental spectral resolution or 
Doppler broadening of the galaxy spectrum. Note that this 
depends on the assumed resolution for each model (can be 
checked by typing > b . R. 

basisf ile= Reference of the models to be used. See table 12. Il for a list 
of the models available with the code. 
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2.3 Running the fitting engine: sfit 



This section details the cah to sfit and some of the options available for the 
fit. This is where most of the work is done, and therefore where most of the 
computation time will be spent. 

2.3.1 Inputs 

The first mandatory argument of sfit is the data structure containing the 
spectrum to be fitted and its metadata. This is the result of the command 
spdata=convertjall(' 'myspectrum. fits' ') as explained in Sec. 12.11 The 
second mandatory argument of sfit is the spectral basis of models b, as ob- 
tained from running bRbasis, explained in Sec. 12.21 The simplest STECKMAP 
fit can then be obtained by issuing: 

> x=sf it (spdata,b) 



ITER 


EVAL 


CPU [s] 




FUNC 


GNORM 


STEPLEN 





1 


1.19e-02 


1, 


. 9995698763700508e+04 


8.3e+02 


O.Oe+00 


100 


123 


6.39e-01 


7, 


, 3005684797028195e+01 


5.9e+01 


l.Oe+00 


150 


543 


2.71e+00 


7, 


, 2729881396562234e+01 


1.6e+01 


5.7e-18 


warning ; 


: too many 


function 


evaluations (543) 






ITER 


EVAL 


CPU [s] 




FUNC 


GNORM 


STEPLEN 





1 


4.78e-03 


5, 


. 3944698728148160e+02 


4.3e+05 


O.Oe+00 


82 


462 


2.32e+00 


7, 


, 2666269919203131e+01 


1.6e+01 


8.5e-21 


ITER 


EVAL 


CPU [s] 




FUNC 


GNORM 


STEPLEN 





1 


5.36e-03 


7, 


, 2666269919203131e+01 


1.6e+01 


O.Oe+00 


100 


473 


2.35e+00 


7, 


, 2600584917043065e+01 


3.4e+01 


l.Oe+00 


111 


507 


2.54e+00 


7, 


, 2599790146012424e+01 


6.5e+01 


1.3e-04 


warning ; 


: too many 


function 


evaluations (507) 








3 14001 10 




5 







1: array (double, 7507) min=0. 783479 max=0. 845731 avg=0. 817507 std=0. 01788 
ebv_star = 0.00900828 



> 

A typical output of sfit is shown above, allowing the user to monitor the fitting 
process. The values of the objective function are shown in the middle column, 
as the minimisation proceeds. This is usuallty close to the of the fit. In this 
example, the fit is bad (the final = 72.6), and the user would attempt to 
improve it. Guidelines for this are given in Sec. SI Turning on the kinematics 
with the option kin=l usually makes a huge difference. 
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2.3.2 Options 



A number of options arc available to enhance the fit, the most important one 
being the kinematics switch (kin=l). All the options can be seen by typing: 

> info.sfit 

f unc sf it (gallist , base , &_ki , &_resi , nMC= , MC= , verb= , maxit= , maxitMC= , noskip= , 
meval= , mevalMC= , meval3= , kin= , epar = , MASK= , sadgues s= , gues s= , mus= , vl im= , 
losguesswidth=,sav=,asc=,RMASK=,zliin=,rmin=,bf =,padv=,f rtoll=,N=,nor=,nde=, 
s igdef = , L 1= , L2= , L3= , L4= , muc= , co= , pl= , muco v= , co v= , parage= , RG= , muv= , mux= , muz= , 

mue= ,mubl= ,mub2= , z2d= , bnpec= , splin= , sigm= , AMRp= , muAMRp= , rrAMRp= , Nb= , lsp= , 

plavgW= , nedgeinask= , dMC= , MC_noise= , pr= , f g= , hidepaths= , piecewise_sad=) 

The most useful options are described in the following table. Others should 
be ignored and not played with. 
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1 to allow for non-parametric LOSVD determination. If 
is used then the spectral models in the basis b must already 
have the correct broadening and shift. 
Velocity range on which the LOSVD is computed. An 
array of 2 elements must be provided, containing the 
lower and upper limits of the velocity range, such as 
vlim= [-1500 . ,1500.], in km/s. The values given as ex- 
amples here are the default. 

if then spectra that already have a results file will not 
be processed. Useful when processing large datasets. 1 
overrides this and reprocesses the spectrum even if it has 
already been processed before. 

type of extinction applied to the model. 1: extinction law 
parametered by E(B-V) taken from [1]. 3: non-parametric 
extinction curve with nde anchor points evenly distributed 
along the spectrum. This is usually a good choice because 
it absorbs smooth flux calibrations errors as well or errors 
in the correction of the instrument responds, default is 3. 
Smoothing parameter for the SAD. (/i^ in the STECMAP 
paper). 

Smoothing parameter for the AMR. (fiz in the STECMAP 
paper). 

Smoothing parameter for the LOSVD. (/iv in the 
STECKMAP paper). 

Smoothing parameter for the non-parametric extinction 
curve. 

Number of anchor points for the non-parametric extinction 
curve. 

1st Binding parameter for the AMR (^c in the STECMAP 
paper). Do not modify. Wildly unstable behaviour may 
ensue. 

2nd Binding parameter for the AMR (^c in the STECMAP 
paper). Do not modify. Wildly unstable behaviour may 
ensue. 

Penalization for the SAD. Allowed values are ' 'Dl' ' (gra- 
dient), "D2" (Laplacian, as in STECMAP paper, or 
' 'D3' 3rd order. Default is Laplacian. 
Penalisation for the LOSVD. Same comments as for LI. 
Penalisation for the AMR. Same comments as above except 
default is gradient. 

Penalisation for the non-parametric extinction curve. Same 
comments as above, default is Laplacian. 
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Chapter 3 

Reading the output files 



When the option sav= 1 is given to sfit, the results are saved. The various results 
files are created in the same directory as the initial data file, upon success of the 
fitting procedure. The .pdb file is the result of the convert_all procedure. It is 
a binary file containing the initial data, a wavelength vector, an error spectrum 
and sometimes a mask if it has been supplied in the proper way (such as when 
dealing with SDSS data). This chapter describes the other files available. 

3.1 Description of the results files 

The format of the files is text as follows: 

1. variable name 

2. length of vector 

3. values 

4. variable name 

5. length of vector 

6. ... 

3.1.1 Stellctr content 

These are the files: res-SAD, res-MASS, res-SFR and res- AMR. They look like 
this (here for the res-MASS file): 

Ages (Myr) 10 



10.0000 
682.190 



23.2692 
1587.40 



54.1455 
3693.75 



125.992 
8595.06 



293.173 
20000.0 



Masses in each time 



bin 



10 



39.4829 
5 . 72450 



17.2968 
251.843 



5 . 35236 
2432.88 



2.45517 
16218.1 



0.997133 
39665.3 
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The variables relative to the stellar content all refer to a time axis, whieh 
can be seen as lookback time or the age of the stellar component. It is always 
given first in these files. The (Myr) indicates the time axis is in Myr. Then 
comes the field of interest, here the masses. Here is a short description of the 
content of each file: 

• res-SAD: The Stellar Age Distribution, i.e. the contribution in flux of 
each component to the observed spectrum. It is normalized so that the 

sum of the SAD is 1. It is the basic quantity STECKMAP is working with 
(rather than masses or SFR). In a typical STECKMAP run the SAD is 
shown in the top plot of the results panel. 

• res-AMR: The Age-Metallicity Relation. This gives the metallicity Z{i) 
of each component i of age age{i), hence defining effectively an age- 
metallicity relation. In extragalactic astronomy, the metallicity is usually 
understood as a total abundance of metals, i.e. the fractional mass of 
everything which is not H or He. With this convention, in the models 
used here the metallicity of the sun is 0.02. Some authors will use 0.019 
instead, but in practice this does not make a huge difference in the inter- 
pretation of the results, in particular as long as one does not embark into 
computing abundances of individual elements. In a typical STECKMAP 
run the AMR is shown as the middle plot of the results panel. 

• res-MASS: This file gives the stellar mass as a function of the age. Note 
that it is understood as an "initial mass", i.e. the mass the given compo- 
nent had at the time of its birth. No mass loss is taken into account here. 
The masses for each time bin i are simply computed as: 

, N sad(i) 
mass{i) = — — (3.1) 

^ ^ M/Liage{i), Z{i)) ^ ' 

where M/L{age,Z) is the mass to light ratio of a stellar population of 
age age and metallicity Z. Again since no mass loss is taken into account 
this refers to the initial m,ass, i.e. it only accounts for the dimming of the 
population, not for the decrease/recycling of its stellar mass. The masses 
are given in principle in solar masses but the normalization is arbitrary: 
the relative masses (i.e. between the various bins) are correct but the 
actual normalization results from the fact that sad is normalized to unity, 
i.e. J2i sad{i) = 1. 

• res-SFR: The Star Formation Rate (SFR) as a function of lookback time. 
It is obtained as: 

sm,-^. ,3.) 

where At{i) is the duration or extent of the age bin i. Note that this is 
not t(i + 1) — t(i) but is rather computed as like At{i) = t{i + 1/2) — t(i — 
1/2). The SFR is given in Mq/ut but is not physically normalized. As 
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for the masses, the relative variations are eorreet but not their absohite 
values, and the actual normalization is essentially inherited from that of 
the masses mass. 



3.1.2 Kinematics 

The broadening function of the data results from the convolution of the LOSVD 
and the instrument's PSF or LSF. However, I will often abusively call it the 
LOSVD, and talk about kinematics eventhough I should be talking about kine- 
matics + instrumental characteristics. The broadening function is stored in the 
res-LOSVD.txt file. It begins as follows: 

V (km/s) 369 

-1003.75 -998.323 -992.892 -987.460 -982.029 

-976.597 -971.165 -965.733 -960.301 -954.869 



This describes the velocity range (in km/s) on which the LOSVD has been 
computed. Next come the actual values of the broadening function for each of 
these velocities: 

g(v) 369 

1.57828e-06 1.62290e-06 1.65155e-06 1.65992e-06 1.63955e-06 
1.61256e-06 1.54141e-06 1.47136e-06 1.33613e-06 1.13671e-06 



The LOSVD or BF g{v) is normalized to unity, i.e. its sum is equal to 1. 



3.1.3 Spectra 

A number of spectra are available in the res-spectra.txt file. In order of 
appearance are the wavelength range (in Angstrom), the data, the best fit- 
ting model, the weight vector, and, if the option cpar=3 has been used, the 
non-parametric extinction curve (maybe calling it non-parametrically adjusted 
continuum (i.e. NPAC would make more sense...). Although these are readily 
plotted in a usual STECKMAP run, this results file should allow the user more 
fiexibility in displaying/ examining his fit. 
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Chapter 4 

Making sense of the results 



This chapter is stiU work in progress. It gives a number of general recipes to 
improve the significance of the result, mostly by improving the fit and running 
Monte Carlo tests. 

4.1 How to improve your fit ? 

While a good (i-e. close to 1) does not necessarily mean that the derived 
SFH is the true SFH of the observed object, it is quite obvious that a SFH 
resulting from a bad fit, i.e. with a > 1-2 for instance (as a rule of thumb), 
can not be trusted. Hence, part of the process of using STECKMAP involves 
improving the fit, or lowering the oi the fit. In some cases, the S/N of the 
data has been overestimated by the observer. It can be set manually using the 
SNRO= keyword convert_all (see Sec. 12. In such a case the user should also 
check afterwards that the residuals of the fit do indeed look like white noise. 
The 2 other important things to investigate are the bad pixels or emission lines, 
and the hyperparameters. 

4.1.1 Masking bad pixels/emission lines 

Since the model libraries only contains stellar models with absorption lines, care 
must be taken in masking the nebular emission lines, possible sky subtraction 
and other bad pixels. A list of common emission lines in the optical is given in 
Appendix A. A mask can be easily generated as a vector of pairs of wavelengths. 
The following mask is designed for masking a 10 Awide region around and 
Hq, for a spectrum that in its restframe. The call to sf it using the mask is also 
shown. 

> mask= [ [4856 . , 4866 . ] , [6558 . , 6568 . ] ] 

> x=sf it(a,b,RMASK=mask) 
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4.1.2 Dependence on hyperparameters 

The default smoothing parameters usuaUy achieve a good balance between keep- 
ing the solution smooth while still fitting the data well. From there, increasing 
mux for instance will usually yield a smoother solution, at the price of a slightly 
larger x^, as shown in [71 [5]. On the other hand, lowering mux will improve the 
is a negligible way, while making the SAD very unstable and sensitive to noise 
(as MonteCarlo experiments show). Lowering mux further will at some point 
result in catastrophic failure of the minimisation procedure, and will return a 
totally chaotic SAD, with a bad x^. It may be informative for the user to play 
around with the smoothing parameters (starting with mux) in a range of a few 
decades (10~^ — 10^) around the defaut values. It will also be useful to check 
out which values have been used in existing publications using STECKMAP 
(given in the foreword to this guide) . 

4.1.3 Summary of guidelines for improving the fit 

In order of importance the following steps must be taken. 

1. Kinematics: if the kinematics are off (kin=0) and the redshift and BF 
are incorrect, the fit can not be good. Turn the kinematics on (kin=l) or 
correct manually the redshift and/or the BF (the models can be broadened 
to an arbitrary resolution using the keyword R= in bRbasisS. 

2. Bad pixels/emission lines. 

3. Re-estimation of the S/N of the data: this is usually a delicate task, 
although in some cases it can be obvious that the S/N is overestimated. 

4. Investigate dependence on hyperparameters. Usually the default parame- 
ters are fine but in some very rare cases better can be obtained with 
non-default values. 

4.2 Assessing the significance of the results 

It is natural and wise to wonder about the sgnificance of our results. There are 
at least 3 types of errors we should care about: 

• Noise in the data: photon noise, read-out noise, background noise... 

• Methodological errors: altough authors of spectral inversion methods do 
not like to think about it there can be systematic effects due to the method 
itself and the choices that have been made in terms of interpolating be- 
tween the ssp models, their normalization etc.... Perhaps the best way to 
investigate these is to use different methods in simple, controlled experi- 
ments, as was done in [1]. 



18 



• Theoretical uncertainties: the ssp models are at the core of galaxy spec- 
tra inversion techniques, and they are not perfect (uncertainties in the 
stellar parameters and coverage of stellar libraries, uncertainties in stellar 
tracks...). A good way to investigate the effect of these is to always analyse 
a given spectrum with a collection of models rather than just one model. 
Several popular models are already plugged in the distributed version of 
STECKMAP (see table 12. ip , and can be used simply by changing the pa- 
rameter basisf ile in the bRbasisS call. This approach was also used in 

In the following we will focus on the methods available in STECKMAP to 
tackle the first item and some of the second, via MonteCalro experiments. 

4.2.1 Robustness of the results: running Monte Carlo ex- 
periments 

Here we understand Monte Carlo experiments as: rerunning the inversion on the 
dataset with a randomized element. The randomized element in STECKMAP 
can either be the noise in the data or the first guess (i.e.e the starting point of 
the minimization method). 

4.2.2 Robustness wrt noise in the data 

Here we want to invert mock data, which will be either the original data with a 
different realization of the noise, or the best model we have of that data, with 
noise added according to the noise vector. The relevant options for sf it are 
listed in Tab. O] 

A typical call will look like: 
x=sf it(a,b,kin=l,noskip=l,sav=l,nMC=5,dMC=' 'bestmodel' ' ,MCjioise=' 'yes' ') 

4.2.3 Robustness wrt first guess choice 

All iterative optimization methods require a starting point, which is usually re- 
ferred to as the first guess. In the case of SFH reconstruction from spectra, scien- 
tists aim at designing a method and an implemetation of it, that will be as little 
as possible sensitive to the first guess. Ideally, we would like STECKMAP to give 
always the same answer, no matter what the first guess was. In STECKMAP, 
the default first guess is a flat distribution for all fields. This flat first guess has 
2 good properties. 

• Convergence: the method usually converges to a good solution (in the 

and regularization sense), which is NOT TRUE for all first guesses: 
in fact, it is fairly easy to come up with a first guess that leads to a 
secondary minimum with a very bad x^, or which leads to divergence of 
the minimization method. 

• Uniformity 
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nMC= 


number of MC experiments to be performed 


dMC= 


determines the mock data to be used for the MC experiments. Two 
choices are possible only: 

dMC=' 'bestmodel' ' will create a mock spectrum by taking the best 
model (corresponding to the derived solution SAD, AMR, BF etc..) 
and noising it using the noise spectrum. The latter has been provided 
by the user during the conversion with convert_all. See Sec. 12.11 for 
the specifications of the noise spectrum or global signal to noise ratio. 
dMC= ' ' data' ' . In this case the mock spectrum will be the original data 
with noise added according to the noise spectrum or the global SNR. 
Note that with dMC=' 'bestmodel' ' we expect ~ 1 for the MC runs, 
while in the case dMC=' 'data' ' , we expect « 2 because the original 
data already contains noise. 


MC_noise= 
MC_noise= 


' ' yes ' ' : mock data is noised according to noise spectrum or global SNR 
' ' no ' ' : mock data is not noised (if mock data is the original data then 
no extra noise is added, same for best model) 



Table 4.1: Options for sf it for performing MonteCarlo experiments to test the 
robustness of the results with respect to noise in the data. 



From my own experience using STECKMAP I acquired confidence in the fact 
that the flat first guess usually does a good job. However there is no proof that 
it is always the best or that it will not sometimes lead to misfits or secondary 
minima. Note that this is true for all available methods OH]. To cope with 
this, many authors note the fact that in order to make sure that they reach a 
"global" minimum and not a secondary minimum, they have to run the fitting 
engine several times, using a family of first guesses. 

Several techniques to vary the first guess have been implemented in sf it. 
They are activated and controlled using the options of Tab. 14.21 They complete 
the options given in Tab. 14.11 which affect the mock data. The preferred setup 
for first guess variation experiments should be : dMC= ' ' data ' ' , MC_noise= ' ' no ' ' ) . 

A typical call in order to test 5 random bursts as first guesses, fitting the 
original data (and not some mock data), will look like: 

x=sfit(a,b,kin=l,noskip=l,sav=l,nMC=5,dMC=' 'data' ' ,fg=' 'RB" ,MCjioise=' 'no") 

4.2.4 Output of MonteCarlo experiments 

The results of each experiments are conveniently stored in the -res . txt files, 
using the same format as the non-MC results, so that for a nMC=5 run, one 
has 6 SADs, AMRs, etc... At the moment the corresponding first guesses are 
not stored in a txt file (altho thats the way it should be!). Here are a few 
lines of example showing how to recover them. In this example all the fields 
are concatenated and not properly normalized. It is nonetheless enough to see 
where the first guess SADs peaks. 

#include "STECKMAP/Pierre/POP/sf it . i" 
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nMC= 


number of MC experiments to be performed 


fg= 




determines the type of first guess to be used for tlie MC experiments. 
Three choices are possible: 


fg=' 


'flat' ' 


will create a flat first guess. This is the default. 


fg=> 


'RB" 


stands for "random burst" . The first guess for the SAD is a single burst 
at a random age within the ages of the SSP basis. The width of the burst 
is 3 age bins. The AMR is flat as well with a random mctallicity within 
the Z range of the SSP basis. The color excess is or the NPEC vector 
is 1 everywhere, depending on the option chosen to deal with extinction, 
and the BF is a gaussian centered on km/s. 


fg=> 


'RG' ' 


stands for "random guess" . All fields are randomized. Not very stable. 
Prefer the 2 first alternatives. 



Table 4.2: Options of sf it for performing MonteCarlo experiments to test the 
robustness of the results with respect to variations of the first guess. 



a=convert_all(fV,zO=0.0021) ; 
b=bRbasis3(,basisf ile="MILES") ; 

x=sfit(a,b,noskip=l,kin=l,nMC=5,dMC="data",fg="RB",MC_noise="no",sav=l,ineval=500) ; 

// plot the first guesses (all fields concatenated) 

upload, a. resfile(l) ; 

ws,l; 

plb.RBs; 
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Chapter 5 

Tweaking STECKMAP 



In some cases you will want to use STECKMAP in a way that was not initially 
planned by the author. This chapter gives a few tricks and corresponding yorick 
codes to do so. 

5.1 Adding a constraint on the AMR 

5.2 Fitting the spectrum with SSPs + power- 
law 
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Appendix A 

Common emission lines in 
the optical 

Emission lines of a variety of chemical elements are usually seen in optical 
spectra of galaxies. Hcire is a compilation of the most common ones. It is 
not meant to be exhaustive but will help the user to design masks efficiently. 
Wavelengths arc given in the air. 



25 



Emission Wavelength Ref. 



line 


(air, A) 




Ha 


6563 


(2) 


H/3 


4861 


(2) 


H7 


4340 


(2) 


R6 


4101 


(2) 


He 


3970 


(2) 


HC 


3889 


(2) 


Hel 


5876 


(1) 


[0 1] 


6300 


(1) 


[0 II] 


3727 


(1) 


[0 III] 


4363 


(1) 


[0 III] 


4959 


(1) 


[0 HI] 


5007 


(1) 


[Nil] 


6584 


(1) 


[Nil] 


6548 


(1) 


[Ne III] 


3869 


(1) 


[S II] 


6716 


(1) 


[S II] 


6731 


(1) 



Table A.l: Common emission lines in optical spectra. Refer- 
ences: (1): http://www.mpa-garching.mpg.de/SDSS/DR7/raw_data.litml, (2): 
http:/ / en.wikipedia.org/ wiki/Balmer .series 
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